
#include "hexahedron.h"
#include "hexahedron20.h"


#include "basis/containers/serializededata.h"
#include "esinfo/eslog.h"
#include "esinfo/meshinfo.h"
#include "mesh/element.h"
#include "math/matrix.dense.h"

using namespace espreso;

void Hexahedron20::setGaussPointsForOrder(int order)
{
	std::vector<double> r, s, t, w;
	if (!Hexahedron::gpw(order, r, s, t, w)) {
		eslog::internalFailure("cannot set Gauss points for a given order.\n");
	}

	N->clear(); N->resize(s.size(), MatrixDense(1, nodes));
	dN->clear(); dN->resize(s.size(), MatrixDense(1, nodes));
	weighFactor->assign(w.begin(), w.end());

	for (size_t i = 0; i < r.size(); i++) {
		MatrixDense &m = (*N)[i];

		// basis function
		m(0, 0) = 0.125 * ((1.0 - r[i]) * (1.0 - s[i]) * (1.0 - t[i]) * (-r[i] - s[i] - t[i] - 2.0));
		m(0, 1) = 0.125 * ((1.0 + r[i]) * (1.0 - s[i]) * (1.0 - t[i]) * ( r[i] - s[i] - t[i] - 2.0));
		m(0, 2) = 0.125 * ((1.0 + r[i]) * (1.0 + s[i]) * (1.0 - t[i]) * ( r[i] + s[i] - t[i] - 2.0));
		m(0, 3) = 0.125 * ((1.0 - r[i]) * (1.0 + s[i]) * (1.0 - t[i]) * (-r[i] + s[i] - t[i] - 2.0));
		m(0, 4) = 0.125 * ((1.0 - r[i]) * (1.0 - s[i]) * (1.0 + t[i]) * (-r[i] - s[i] + t[i] - 2.0));
		m(0, 5) = 0.125 * ((1.0 + r[i]) * (1.0 - s[i]) * (1.0 + t[i]) * ( r[i] - s[i] + t[i] - 2.0));
		m(0, 6) = 0.125 * ((1.0 + r[i]) * (1.0 + s[i]) * (1.0 + t[i]) * ( r[i] + s[i] + t[i] - 2.0));
		m(0, 7) = 0.125 * ((1.0 - r[i]) * (1.0 + s[i]) * (1.0 + t[i]) * (-r[i] + s[i] + t[i] - 2.0));

		m(0, 8) =  0.25 * ((1.0 - pow(r[i], 2.0)) * (1.0 - s[i]) * (1.0 - t[i]));
		m(0, 9) =  0.25 * ((1.0 + r[i]) * (1.0 - pow(s[i], 2.0)) * (1.0 - t[i]));
		m(0, 10) = 0.25 * ((1.0 - pow(r[i], 2.0)) * (1.0 + s[i]) * (1.0 - t[i]));
		m(0, 11) = 0.25 * ((1.0 - r[i]) * (1.0 - pow(s[i], 2.0)) * (1.0 - t[i]));
		m(0, 12) = 0.25 * ((1.0 - pow(r[i], 2.0)) * (1.0 - s[i]) * (1.0 + t[i]));
		m(0, 13) = 0.25 * ((1.0 + r[i]) * (1.0 - pow(s[i], 2.0)) * (1.0 + t[i]));
		m(0, 14) = 0.25 * ((1.0 - pow(r[i], 2.0)) * (1.0 + s[i]) * (1.0 + t[i]));
		m(0, 15) = 0.25 * ((1.0 - r[i]) * (1.0 - pow(s[i], 2.0)) * (1.0 + t[i]));
		m(0, 16) = 0.25 * ((1.0 - r[i]) * (1.0 - s[i]) * (1.0 - pow(t[i], 2.0)));
		m(0, 17) = 0.25 * ((1.0 + r[i]) * (1.0 - s[i]) * (1.0 - pow(t[i], 2.0)));
		m(0, 18) = 0.25 * ((1.0 + r[i]) * (1.0 + s[i]) * (1.0 - pow(t[i], 2.0)));
		m(0, 19) = 0.25 * ((1.0 - r[i]) * (1.0 + s[i]) * (1.0 - pow(t[i], 2.0)));
	}

	for (size_t i = 0; i < r.size(); i++) {
		MatrixDense &m = (*dN)[i];

		m(0, 0) =  ((s[i] - 1.0) * (t[i] - 1.0) * (r[i] + s[i] + t[i] + 2.0)) / 8.0 + ((r[i] - 1.0) * (s[i] - 1.0) *         (t[i] - 1.0)) / 8.0;
		m(0, 1) =  ((r[i] + 1.0) * (s[i] - 1.0) *         (t[i] - 1.0)) / 8.0 - ((s[i] - 1.0) * (t[i] - 1.0) * (s[i] - r[i] + t[i] + 2.0)) / 8.0;
		m(0, 2) = -((s[i] + 1.0) * (t[i] - 1.0) * (r[i] + s[i] - t[i] - 2.0)) / 8.0 - ((r[i] + 1.0) * (s[i] + 1.0) *         (t[i] - 1.0)) / 8.0;
		m(0, 3) = -((s[i] + 1.0) * (t[i] - 1.0) * (r[i] - s[i] + t[i] + 2.0)) / 8.0 - ((r[i] - 1.0) * (s[i] + 1.0) *         (t[i] - 1.0)) / 8.0;
		m(0, 4) = -((s[i] - 1.0) * (t[i] + 1.0) * (r[i] + s[i] - t[i] + 2.0)) / 8.0 - ((r[i] - 1.0) * (s[i] - 1.0) *         (t[i] + 1.0)) / 8.0;
		m(0, 5) = -((s[i] - 1.0) * (t[i] + 1.0) * (r[i] - s[i] + t[i] - 2.0)) / 8.0 - ((r[i] + 1.0) * (s[i] - 1.0) *         (t[i] + 1.0)) / 8.0;
		m(0, 6) =  ((s[i] + 1.0) * (t[i] + 1.0) * (r[i] + s[i] + t[i] - 2.0)) / 8.0 + ((r[i] + 1.0) * (s[i] + 1.0) *         (t[i] + 1.0)) / 8.0;
		m(0, 7) =  ((s[i] + 1.0) * (t[i] + 1.0) * (r[i] - s[i] - t[i] + 2.0)) / 8.0 + ((r[i] - 1.0) * (s[i] + 1.0) *         (t[i] + 1.0)) / 8.0;

		m(0, 8)  = -(r[i] * (s[i] - 1.0) * (t[i] - 1.0)) / 2.0;
		m(0, 9)  =  ((pow(s[i], 2.0) - 1.0) * (t[i] - 1.0)) / 4.0;
		m(0, 10) =  (r[i] * (s[i] + 1.0) * (t[i] - 1.0)) / 2.0;
		m(0, 11) = -((pow(s[i], 2.0) - 1.0) * (t[i] - 1.0)) / 4.0;
		m(0, 12) =  (r[i] * (s[i] - 1.0) * (t[i] + 1.0)) / 2.0;
		m(0, 13) = -((pow(s[i], 2.0) - 1.0) * (t[i] + 1.0)) / 4.0;
		m(0, 14) = -(r[i] * (s[i] + 1.0) * (t[i] + 1.0)) / 2.0;
		m(0, 15) =  ((pow(s[i], 2.0) - 1.0) * (t[i] + 1.0)) / 4.0;
		m(0, 16) = -((pow(t[i], 2.0) - 1.0) * (s[i] - 1.0)) / 4.0;
		m(0, 17) =  ((pow(t[i], 2.0) - 1.0) * (s[i] - 1.0)) / 4.0;
		m(0, 18) = -((pow(t[i], 2.0) - 1.0) * (s[i] + 1.0)) / 4.0;
		m(0, 19) =  ((pow(t[i], 2.0) - 1.0) * (s[i] + 1.0)) / 4.0;


		m(1, 0) =  ((r[i] - 1.0) * (t[i] - 1.0) * (r[i] + s[i] + t[i] + 2.0)) / 8.0 + ((r[i] - 1.0) * (s[i] - 1.0) *         (t[i] - 1.0)) / 8.0;
		m(1, 1) = -((r[i] + 1.0) * (t[i] - 1.0) * (s[i] - r[i] + t[i] + 2.0)) / 8.0 - ((r[i] + 1.0) * (s[i] - 1.0) *         (t[i] - 1.0)) / 8.0;
		m(1, 2) = -((r[i] + 1.0) * (t[i] - 1.0) * (r[i] + s[i] - t[i] - 2.0)) / 8.0 - ((r[i] + 1.0) * (s[i] + 1.0) *         (t[i] - 1.0)) / 8.0;
		m(1, 3) =  ((r[i] - 1.0) * (s[i] + 1.0) *         (t[i] - 1.0)) / 8.0 - ((r[i] - 1.0) * (t[i] - 1.0) * (r[i] - s[i] + t[i] + 2.0)) / 8.0;
		m(1, 4) = -((r[i] - 1.0) * (t[i] + 1.0) * (r[i] + s[i] - t[i] + 2.0)) / 8.0 - ((r[i] - 1.0) * (s[i] - 1.0) *         (t[i] + 1.0)) / 8.0;
		m(1, 5) =  ((r[i] + 1.0) * (s[i] - 1.0) *         (t[i] + 1.0)) / 8.0 - ((r[i] + 1.0) * (t[i] + 1.0) * (r[i] - s[i] + t[i] - 2.0)) / 8.0;
		m(1, 6) =  ((r[i] + 1.0) * (t[i] + 1.0) * (r[i] + s[i] + t[i] - 2.0)) / 8.0 + ((r[i] + 1.0) * (s[i] + 1.0) *         (t[i] + 1.0)) / 8.0;
		m(1, 7) =  ((r[i] - 1.0) * (t[i] + 1.0) * (r[i] - s[i] - t[i] + 2.0)) / 8.0 - ((r[i] - 1.0) * (s[i] + 1.0) *         (t[i] + 1.0)) / 8.0;

		m(1, 8)  = -((pow(r[i], 2.0) - 1.0) * (t[i] - 1.0)) / 4.0;
		m(1, 9)  =  (s[i] * (r[i] + 1.0) * (t[i] - 1.0)) / 2.0;
		m(1, 10) =  ((pow(r[i], 2.0) - 1.0) * (t[i] - 1.0)) / 4.0;
		m(1, 11) = -(s[i] * (r[i] - 1.0) * (t[i] - 1.0)) / 2.0;
		m(1, 12) =  ((pow(r[i], 2.0) - 1.0) * (t[i] + 1.0)) / 4.0;
		m(1, 13) = -(s[i] * (r[i] + 1.0) * (t[i] + 1.0)) / 2.0;
		m(1, 14) = -((pow(r[i], 2.0) - 1.0) * (t[i] + 1.0)) / 4.0;
		m(1, 15) =  (s[i] * (r[i] - 1.0) * (t[i] + 1.0)) / 2.0;
		m(1, 16) = -((pow(t[i], 2.0) - 1.0) * (r[i] - 1.0)) / 4.0;
		m(1, 17) =  ((pow(t[i], 2.0) - 1.0) * (r[i] + 1.0)) / 4.0;
		m(1, 18) = -((pow(t[i], 2.0) - 1.0) * (r[i] + 1.0)) / 4.0;
		m(1, 19) =  ((pow(t[i], 2.0) - 1.0) * (r[i] - 1.0)) / 4.0;

		m(2, 0) =  ((r[i] - 1.0) * (s[i] - 1.0) * (r[i] + s[i] + t[i] + 2.0)) / 8.0 + ((r[i] - 1.0) * (s[i] - 1.0) *         (t[i] - 1.0)) / 8.0;
		m(2, 1) = -((r[i] + 1.0) * (s[i] - 1.0) * (s[i] - r[i] + t[i] + 2.0)) / 8.0 - ((r[i] + 1.0) * (s[i] - 1.0) *         (t[i] - 1.0)) / 8.0;
		m(2, 2) =  ((r[i] + 1.0) * (s[i] + 1.0) *         (t[i] - 1.0)) / 8.0 - ((r[i] + 1.0) * (s[i] + 1.0) * (r[i] + s[i] - t[i] - 2.0)) / 8.0;
		m(2, 3) = -((r[i] - 1.0) * (s[i] + 1.0) * (r[i] - s[i] + t[i] + 2.0)) / 8.0 - ((r[i] - 1.0) * (s[i] + 1.0) *         (t[i] - 1.0)) / 8.0;
		m(2, 4) =  ((r[i] - 1.0) * (s[i] - 1.0) *         (t[i] + 1.0)) / 8.0 - ((r[i] - 1.0) * (s[i] - 1.0) * (r[i] + s[i] - t[i] + 2.0)) / 8.0;
		m(2, 5) = -((r[i] + 1.0) * (s[i] - 1.0) * (r[i] - s[i] + t[i] - 2.0)) / 8.0 - ((r[i] + 1.0) * (s[i] - 1.0) *         (t[i] + 1.0)) / 8.0;
		m(2, 6) =  ((r[i] + 1.0) * (s[i] + 1.0) * (r[i] + s[i] + t[i] - 2.0)) / 8.0 + ((r[i] + 1.0) * (s[i] + 1.0) *         (t[i] + 1.0)) / 8.0;
		m(2, 7) =  ((r[i] - 1.0) * (s[i] + 1.0) * (r[i] - s[i] - t[i] + 2.0)) / 8.0 - ((r[i] - 1.0) * (s[i] + 1.0) *         (t[i] + 1.0)) / 8.0;

		m(2, 8)  = -((pow(r[i], 2.0) - 1.0) * (s[i] - 1.0)) / 4.0;
		m(2, 9)  =  ((pow(s[i], 2.0) - 1.0) * (r[i] + 1.0)) / 4.0;
		m(2, 10) =  ((pow(r[i], 2.0) - 1.0) * (s[i] + 1.0)) / 4.0;
		m(2, 11) = -((pow(s[i], 2.0) - 1.0) * (r[i] - 1.0)) / 4.0;
		m(2, 12) =  ((pow(r[i], 2.0) - 1.0) * (s[i] - 1.0)) / 4.0;
		m(2, 13) = -((pow(s[i], 2.0) - 1.0) * (r[i] + 1.0)) / 4.0;
		m(2, 14) = -((pow(r[i], 2.0) - 1.0) * (s[i] + 1.0)) / 4.0;
		m(2, 15) =  ((pow(s[i], 2.0) - 1.0) * (r[i] - 1.0)) / 4.0;
		m(2, 16) = -(t[i] * (r[i] - 1.0) * (s[i] - 1.0)) / 2.0;
		m(2, 17) =  (t[i] * (r[i] + 1.0) * (s[i] - 1.0)) / 2.0;
		m(2, 18) = -(t[i] * (r[i] + 1.0) * (s[i] + 1.0)) / 2.0;
		m(2, 19) =  (t[i] * (r[i] - 1.0) * (s[i] + 1.0)) / 2.0;
	}
}

void Hexahedron20::setBaseFunctions(Element &self)
{
	size_t GPCount = 8, nodeCount = 20;

	std::vector<std::vector<double> > rst(3);

	switch (GPCount) {
	case 8: {
		double v = 0.577350269189625953;
		rst[0] = {  v,  v,  v,  v, -v, -v, -v, -v };
		rst[1] = { -v, -v,  v,  v, -v, -v,  v,  v };
		rst[2] = { -v,  v, -v,  v, -v,  v, -v,  v };
		break;
	}
	case 14: {
		double v1 = 0.758786910639329015;
		double v2 = 0.795822425754222018;
		double v3 = 0;
		rst[0] = { -v1,  v1,  v1, -v1, -v1,  v1,  v1, -v1,  v3,  v3,  v2, v3, -v2, v3 };
		rst[1] = { -v1, -v1,  v1,  v1, -v1, -v1,  v1,  v1,  v3, -v2,  v3, v2,  v3, v3 };
		rst[2] = { -v1, -v1, -v1, -v1,  v1,  v1,  v1,  v1, -v2,  v3,  v3, v3,  v3, v2 };
		break;
	}
	default:
		exit(1);
	}

	self.N = new std::vector<MatrixDense>(GPCount, MatrixDense(1, nodeCount));
	self.dN = new std::vector<MatrixDense>(GPCount, MatrixDense(3, nodeCount));
	self.weighFactor = new std::vector<double>();

	for (unsigned int i = 0; i < GPCount; i++) {
		MatrixDense &m = (*self.N)[i];

		double r = rst[0][i];
		double s = rst[1][i];
		double t = rst[2][i];

		// basis function
		m(0, 0) = 0.125 * ((1.0 - r) * (1.0 - s) * (1.0 - t) * (-r - s - t - 2.0));
		m(0, 1) = 0.125 * ((1.0 + r) * (1.0 - s) * (1.0 - t) * ( r - s - t - 2.0));
		m(0, 2) = 0.125 * ((1.0 + r) * (1.0 + s) * (1.0 - t) * ( r + s - t - 2.0));
		m(0, 3) = 0.125 * ((1.0 - r) * (1.0 + s) * (1.0 - t) * (-r + s - t - 2.0));
		m(0, 4) = 0.125 * ((1.0 - r) * (1.0 - s) * (1.0 + t) * (-r - s + t - 2.0));
		m(0, 5) = 0.125 * ((1.0 + r) * (1.0 - s) * (1.0 + t) * ( r - s + t - 2.0));
		m(0, 6) = 0.125 * ((1.0 + r) * (1.0 + s) * (1.0 + t) * ( r + s + t - 2.0));
		m(0, 7) = 0.125 * ((1.0 - r) * (1.0 + s) * (1.0 + t) * (-r + s + t - 2.0));

		m(0, 8) =  0.25 * ((1.0 - pow(r, 2.0)) * (1.0 - s) * (1.0 - t));
		m(0, 9) =  0.25 * ((1.0 + r) * (1.0 - pow(s, 2.0)) * (1.0 - t));
		m(0, 10) = 0.25 * ((1.0 - pow(r, 2.0)) * (1.0 + s) * (1.0 - t));
		m(0, 11) = 0.25 * ((1.0 - r) * (1.0 - pow(s, 2.0)) * (1.0 - t));
		m(0, 12) = 0.25 * ((1.0 - pow(r, 2.0)) * (1.0 - s) * (1.0 + t));
		m(0, 13) = 0.25 * ((1.0 + r) * (1.0 - pow(s, 2.0)) * (1.0 + t));
		m(0, 14) = 0.25 * ((1.0 - pow(r, 2.0)) * (1.0 + s) * (1.0 + t));
		m(0, 15) = 0.25 * ((1.0 - r) * (1.0 - pow(s, 2.0)) * (1.0 + t));
		m(0, 16) = 0.25 * ((1.0 - r) * (1.0 - s) * (1.0 - pow(t, 2.0)));
		m(0, 17) = 0.25 * ((1.0 + r) * (1.0 - s) * (1.0 - pow(t, 2.0)));
		m(0, 18) = 0.25 * ((1.0 + r) * (1.0 + s) * (1.0 - pow(t, 2.0)));
		m(0, 19) = 0.25 * ((1.0 - r) * (1.0 + s) * (1.0 - pow(t, 2.0)));
	}

	for (unsigned int i = 0; i < GPCount; i++) {
		MatrixDense &m = (*self.dN)[i];

		double r = rst[0][i];
		double s = rst[1][i];
		double t = rst[2][i];


		// dNr - derivation of basis function
		m(0, 0) =  ((s - 1.0) * (t - 1.0) * (r + s + t + 2.0)) / 8.0 + ((r - 1.0) * (s - 1.0) *         (t - 1.0)) / 8.0;
		m(0, 1) =  ((r + 1.0) * (s - 1.0) *         (t - 1.0)) / 8.0 - ((s - 1.0) * (t - 1.0) * (s - r + t + 2.0)) / 8.0;
		m(0, 2) = -((s + 1.0) * (t - 1.0) * (r + s - t - 2.0)) / 8.0 - ((r + 1.0) * (s + 1.0) *         (t - 1.0)) / 8.0;
		m(0, 3) = -((s + 1.0) * (t - 1.0) * (r - s + t + 2.0)) / 8.0 - ((r - 1.0) * (s + 1.0) *         (t - 1.0)) / 8.0;
		m(0, 4) = -((s - 1.0) * (t + 1.0) * (r + s - t + 2.0)) / 8.0 - ((r - 1.0) * (s - 1.0) *         (t + 1.0)) / 8.0;
		m(0, 5) = -((s - 1.0) * (t + 1.0) * (r - s + t - 2.0)) / 8.0 - ((r + 1.0) * (s - 1.0) *         (t + 1.0)) / 8.0;
		m(0, 6) =  ((s + 1.0) * (t + 1.0) * (r + s + t - 2.0)) / 8.0 + ((r + 1.0) * (s + 1.0) *         (t + 1.0)) / 8.0;
		m(0, 7) =  ((s + 1.0) * (t + 1.0) * (r - s - t + 2.0)) / 8.0 + ((r - 1.0) * (s + 1.0) *         (t + 1.0)) / 8.0;

		m(0, 8)  = -(r * (s - 1.0) * (t - 1.0)) / 2.0;
		m(0, 9)  =  ((pow(s, 2.0) - 1.0) * (t - 1.0)) / 4.0;
		m(0, 10) =  (r * (s + 1.0) * (t - 1.0)) / 2.0;
		m(0, 11) = -((pow(s, 2.0) - 1.0) * (t - 1.0)) / 4.0;
		m(0, 12) =  (r * (s - 1.0) * (t + 1.0)) / 2.0;
		m(0, 13) = -((pow(s, 2.0) - 1.0) * (t + 1.0)) / 4.0;
		m(0, 14) = -(r * (s + 1.0) * (t + 1.0)) / 2.0;
		m(0, 15) =  ((pow(s, 2.0) - 1.0) * (t + 1.0)) / 4.0;
		m(0, 16) = -((pow(t, 2.0) - 1.0) * (s - 1.0)) / 4.0;
		m(0, 17) =  ((pow(t, 2.0) - 1.0) * (s - 1.0)) / 4.0;
		m(0, 18) = -((pow(t, 2.0) - 1.0) * (s + 1.0)) / 4.0;
		m(0, 19) =  ((pow(t, 2.0) - 1.0) * (s + 1.0)) / 4.0;


		// dNs - derivation of basis function
		m(1, 0) =  ((r - 1.0) * (t - 1.0) * (r + s + t + 2.0)) / 8.0 + ((r - 1.0) * (s - 1.0) *         (t - 1.0)) / 8.0;
		m(1, 1) = -((r + 1.0) * (t - 1.0) * (s - r + t + 2.0)) / 8.0 - ((r + 1.0) * (s - 1.0) *         (t - 1.0)) / 8.0;
		m(1, 2) = -((r + 1.0) * (t - 1.0) * (r + s - t - 2.0)) / 8.0 - ((r + 1.0) * (s + 1.0) *         (t - 1.0)) / 8.0;
		m(1, 3) =  ((r - 1.0) * (s + 1.0) *         (t - 1.0)) / 8.0 - ((r - 1.0) * (t - 1.0) * (r - s + t + 2.0)) / 8.0;
		m(1, 4) = -((r - 1.0) * (t + 1.0) * (r + s - t + 2.0)) / 8.0 - ((r - 1.0) * (s - 1.0) *         (t + 1.0)) / 8.0;
		m(1, 5) =  ((r + 1.0) * (s - 1.0) *         (t + 1.0)) / 8.0 - ((r + 1.0) * (t + 1.0) * (r - s + t - 2.0)) / 8.0;
		m(1, 6) =  ((r + 1.0) * (t + 1.0) * (r + s + t - 2.0)) / 8.0 + ((r + 1.0) * (s + 1.0) *         (t + 1.0)) / 8.0;
		m(1, 7) =  ((r - 1.0) * (t + 1.0) * (r - s - t + 2.0)) / 8.0 - ((r - 1.0) * (s + 1.0) *         (t + 1.0)) / 8.0;

		m(1, 8)  = -((pow(r, 2.0) - 1.0) * (t - 1.0)) / 4.0;
		m(1, 9)  =  (s * (r + 1.0) * (t - 1.0)) / 2.0;
		m(1, 10) =  ((pow(r, 2.0) - 1.0) * (t - 1.0)) / 4.0;
		m(1, 11) = -(s * (r - 1.0) * (t - 1.0)) / 2.0;
		m(1, 12) =  ((pow(r, 2.0) - 1.0) * (t + 1.0)) / 4.0;
		m(1, 13) = -(s * (r + 1.0) * (t + 1.0)) / 2.0;
		m(1, 14) = -((pow(r, 2.0) - 1.0) * (t + 1.0)) / 4.0;
		m(1, 15) =  (s * (r - 1.0) * (t + 1.0)) / 2.0;
		m(1, 16) = -((pow(t, 2.0) - 1.0) * (r - 1.0)) / 4.0;
		m(1, 17) =  ((pow(t, 2.0) - 1.0) * (r + 1.0)) / 4.0;
		m(1, 18) = -((pow(t, 2.0) - 1.0) * (r + 1.0)) / 4.0;
		m(1, 19) =  ((pow(t, 2.0) - 1.0) * (r - 1.0)) / 4.0;

		// dNt - derivation of basis function
		m(2, 0) =  ((r - 1.0) * (s - 1.0) * (r + s + t + 2.0)) / 8.0 + ((r - 1.0) * (s - 1.0) *         (t - 1.0)) / 8.0;
		m(2, 1) = -((r + 1.0) * (s - 1.0) * (s - r + t + 2.0)) / 8.0 - ((r + 1.0) * (s - 1.0) *         (t - 1.0)) / 8.0;
		m(2, 2) =  ((r + 1.0) * (s + 1.0) *         (t - 1.0)) / 8.0 - ((r + 1.0) * (s + 1.0) * (r + s - t - 2.0)) / 8.0;
		m(2, 3) = -((r - 1.0) * (s + 1.0) * (r - s + t + 2.0)) / 8.0 - ((r - 1.0) * (s + 1.0) *         (t - 1.0)) / 8.0;
		m(2, 4) =  ((r - 1.0) * (s - 1.0) *         (t + 1.0)) / 8.0 - ((r - 1.0) * (s - 1.0) * (r + s - t + 2.0)) / 8.0;
		m(2, 5) = -((r + 1.0) * (s - 1.0) * (r - s + t - 2.0)) / 8.0 - ((r + 1.0) * (s - 1.0) *         (t + 1.0)) / 8.0;
		m(2, 6) =  ((r + 1.0) * (s + 1.0) * (r + s + t - 2.0)) / 8.0 + ((r + 1.0) * (s + 1.0) *         (t + 1.0)) / 8.0;
		m(2, 7) =  ((r - 1.0) * (s + 1.0) * (r - s - t + 2.0)) / 8.0 - ((r - 1.0) * (s + 1.0) *         (t + 1.0)) / 8.0;

		m(2, 8)  = -((pow(r, 2.0) - 1.0) * (s - 1.0)) / 4.0;
		m(2, 9)  =  ((pow(s, 2.0) - 1.0) * (r + 1.0)) / 4.0;
		m(2, 10) =  ((pow(r, 2.0) - 1.0) * (s + 1.0)) / 4.0;
		m(2, 11) = -((pow(s, 2.0) - 1.0) * (r - 1.0)) / 4.0;
		m(2, 12) =  ((pow(r, 2.0) - 1.0) * (s - 1.0)) / 4.0;
		m(2, 13) = -((pow(s, 2.0) - 1.0) * (r + 1.0)) / 4.0;
		m(2, 14) = -((pow(r, 2.0) - 1.0) * (s + 1.0)) / 4.0;
		m(2, 15) =  ((pow(s, 2.0) - 1.0) * (r - 1.0)) / 4.0;
		m(2, 16) = -(t * (r - 1.0) * (s - 1.0)) / 2.0;
		m(2, 17) =  (t * (r + 1.0) * (s - 1.0)) / 2.0;
		m(2, 18) = -(t * (r + 1.0) * (s + 1.0)) / 2.0;
		m(2, 19) =  (t * (r - 1.0) * (s + 1.0)) / 2.0;
	}

	switch (GPCount) {
	case 8: {
		self.weighFactor->resize(8, 1.0);
		break;
	}
	case 14: {
		self.weighFactor->resize(8, 0.335180055401662);
		self.weighFactor->resize(14, 0.886426592797784);
		break;
	}
	default:
		exit(1);
	}

	BaseFunctions::created(self);
}


